function Tr = TR(Td, T0, t)
% 式(15)
Tp = (2/3)*Td;
if t > Tp
    Tr = Td;
else
    Tr = Td - (Td - T0)*(t - Tp)^2/(Tp^2);
end